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ABSTRACT 

Robust fast methods to classify variable light curves in large sky surveys 
are becoming increasingly important. While it is relatively straightforward to 
identify common periodic stars and particular transient events (supernovae, no- 
vae, microlensing events), there is no equivalent for non-periodic continuously 
varying sources (quasars, aperiodic stellar variability). In this paper we present 
a fast method for modeling and classifying such sources. We demonstrate the 
method using ~ 86, 000 variable sources from the OGLE-II survey of the LMC 
and ~ 2, 700 mid-IR selected quasar candidates from the OGLE-III survey of 
the LMC and SMC. We discuss the location of common variability classes in the 
parameter space of the model. In particular, we show that quasars occupy a 
distinct region of variability space, providing a simple quantitative approach to 
the variability selection of quasars. 
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Introduction 



The focus of astronomy on the transient universe is steadily increasing because of the 

wealth of astrophysical information provided by time variability. Examp les of past and cur- 

rent v ariabil ity surveys are the Galactic microlensing survey s (e.g., OGLE.Iudalski. Kubiak fc Szymahski 



Becker et al 



19971 : MOA. ISumi et al.ll2003l: M ACHO. lAlcock et al.ll2000l : ER OS. Beaulieu et a 



20051 ). s earches for both local fLOS S. iFilippenko et al 



Miknaitis et al 



per MACHO, 

cosmologically distan t (e.g., SDSS. ISako et al.ll2008l ; CFHTLS. ICabanac et al 



2007) supernovae, all sky variability sur veys (e.g., AS AS. iPoimariski 



QUEST, IVivas et all 12004 ; NSVS, IWozniak et al. 



prom pt 7-ray burst monitors (e.g ., BATS E. iMeegan et al. 



2004; Catalina, Drake et al 



1995; Su- 



200l|) and 



2007; Es sence, 



20021 ; ROTSE. iRvkoff et al.ll2005l : SWIFT. [Burrows et al.ll2005h . A complete, or even partial 



1992; RAPTOR. 



1997 



2009) and 



Vestrand et al. 



listing of more directed monitoring projects would be impossible, and the scale of variability 
surveys continues to grow rapi dly with projects suc h as OGLE- IV (A. Udalski et al. 2010, 
in preparation), Pan-STARRS JHodapp et alil2004h and LSST Jlvezic et al.ll2008h . 



A common problem for variability surveys is the classificatio n of light curves giyen th e 



enormous range of physical processes leading to variability (e.g., ISterken fc Jascheklll996l ). 
Classification is well developed for simple periodic sources (e.g., Cepheids, RR Lyrae, eclips- 
ing binaries, etc.) and impulsive (e.g., microlensing) or explosive sources (e.g., supernovae, 
novae), but less well developed for sources with non-periodic, continuous or very long time- 
scale variability such as quasars. Most efforts at classification have focused on the broad 



different variable classes using a broad ran ge of methods (e.g 



Ever 


2002; 


Pojmahski 


2002 



Ever &; Blakd 120051 ; ISoszyhski et al.l l2008bl ) . M ore sophisticated methods for cla ssification 



of variable objects incl ude neural networks (e.g., iBelokurov. Evans &; Le Dull2004l ) and self- 



ore 



ganizing maps (e.g., Wyrzykowski fc Belokurovl 20081 ). 



Attempts at classif ying non-per i odic, continuously variable sources have largely focused 
on identifying quasars. ISesar et al.l (120071 ) found that essentially all quasars vary at some 
level on long time scales, with 90% of Sloan Digital Sky Survey (SDSS) quasars having 
rms variability > 0.03 mag on a 6 year time scale. The average variability of quasars 
can be describ ed by a structure function with greater varia bility amplitudes on longer time 
baselines (e.g., Vanden Berk et aDl2004 ; de Vries et al. and this has also been observed 



for smaller numbers of individually m onitored quasars (e.g., ICurti et al.lll985l ; iNelson 



Collin et al. 



1983 



2004a 



2002 



Veron Sz Hawkins 



Dobrzycki et al 



1995 



2005 



or 



Sergeev et al.ll2005r). Thus, variability searches for quasa r s (e.g., 



Dobrzycki et al.l l2003al ; iGeha et al.l 120031 ; iRengstorf et al. 



1996; 



Hawkins 



Sumi et al.ll2005l more recently) generally look for non-periodic 



sources with long time scale variability. The most formal approach to date is that used by 
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Eyer (2002) and Sumi et al. (2005), where the selection criterion was based on the slope of 
the structure function of the light curves. These searches have generally been more successful 
in relatively empty extragalactic regions than in those with high stellar densities. 

A significant problem in these attempts is that there was no simple method to reduce a 
non-periodic, continuously variable light curve to a small set of numbers, such as the periods 



and am plitudes of a Fourier series, which can be used for classification. Recently, iKelly et al. 



( 120091 ) found that quasar light curves could be well modeled as a stochastic process with an 
exponential covariance matrix characterized by an amplitude and a time scale based on a 
small sample (~ 100) of quasars with multi-year light curves. Physically, the model is a 
damped random walk, and it has a broken power law structure function consistent with 
studies of quasar structure functions. Unfortunately, there are few large samples of quasars 
with extensive monitoring data. To our knowledge, t here are the quasars in the SPSS equa - 
torial strip, with rou ghly 60 epochs oyer ~si x years (ISesar et al.l 120071 ; iBramich et al.ll2008l ). 
the QUEST survey (IRengstorf et al.ll2004bl ). whose light curves ar e not publicly available. 



and a ny quasars lying in the microlensi ng survey regions. Recently, 



( 120091 ) used mid-IR color selection (e.g. JStern et al.ll2005l . also see lLacy et al 



Kozlows 



a fc Kochanek 
200j) to idem 



tify a large sample of quasar candidates in the Magellanic Clouds. While they have yet to 
be spectroscopically confirmed, the nature of the selection method means that purity of the 
main sample of ca ndidates should be high. By combin i ng these candidates with the OGLE- 
II and OGLE-III JlJdalski. Kubiak fc Szvmahskil fl997l : llJdalski et alihoosl ) light curves, we 
can verify that quasar light curves are well described by a st ochastic process based on a 
much larger sample of objects than used by IKelly et al.l (120091 ) and then explore using this 
approach to light curve modeling for classifying light curves in general, including periodic 
sources. Obviously, a damped random walk is not an optimal model for periodic sources, 
but there is no technical difficulty in modeling them using the stochastic process. 

This paper is organized as follows. In Section [2] we describe the light curve data and 
the existing catalogs of variable sources we can use for classifying light curves. In Sec- 
tion [3] we descr ibe the mathemat ical model. We do not use the particular forecasting 
methodology of Kelly et aP (2009 ). but a more stati s tically optim al approach based on 
Press. Rvbicki fc Hewittl Jl992ah and iRvbicki fc Press! Jl992l . ll994T ) that we detail in the 



appendix. In Sec tion H] w e examine how well the model works for the known quasars from 
Kelly et al.l (120091 ) and the lKozlowski fc Kochanek! (120091 ) quasar candidates, the distribution 
of quasars in the parameter space of the model and the completeness of variability selection. 
In Section \5\ we examine the problem of contamination by determining the distribution of 
the general population of variable LMC sources in the model parameters. We also examine 
how the stochastic process models treat true periodic variables. In Section [6] we combine 
these results and discuss the problem of variability selecting quasars in dense stellar fields 
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like the LMC. We summarize our results and discuss future uses in Section [7J 



Data 



The data used in this paper were collected with the 1.3 m Warsaw telescope at the 
Las Campanas Observatory, Chile during the second (1996-2000) and third (2001- 2008) 
phases of the OGLE project (jUdalski. Kubiak fc Szymariskilll997l ; lUdalski et al.ll2008l ). The 
majority of the OGLE data were collected in the /-band (~ 370 points per light curve in 
OGLE-III) with a smaller number of points in the V^-band (~ 40 points per light curve). 
The OGLE-III light curves for many common variable stars are available on the interneo 
The OGLE- II database of l i ght curves (/-band only) can be ac cessed through several web 



interfaced! (jSzymanskilfeOOa ; lUdalski. Kubiak fc Szvma risk 



OGLE samples of classical C epheids (Cep 



1997). We make use of published 



Soszyhski et al.l 1 2008b), double m ode Cepheids 



(dCep; ISoszyhski et al.ll2008al ). RR Lyrae s tars (RRLvr; ISoszyhski et al.ll2009D . OGLE Small 
Amplitude Variable Red Giants (OSARGs; Soszyhski et al. 2004a), eclipsing binaries^] (ECL; 
Wyrzykowski et al. 2003 ). ellipsoi dal variable red g iant (ELL; ISoszyhski et al. 2004b ). long 



secondary period variables (LSPs 



Soszvhskil 120071). an d long period variables (Miras, LPVs, 



and other semiregular variables, ISoszyhski et al.ll2005f) in the LMC. We separately extracted 
the OGLE-II light curves of Be stars from lKeller et all (booi ). ~ 300 OGLE-II and ~ 2700 



OGLE -III light curves of the mid-IR selected quas ar candidates fromlKozlowski fc Kochanek 
( 120091 ) . We also reanalyzed the 109 quasars from iKelly et al.l (120091 ) . 



In addition to these specific samples we downloaded ~ 150,000 "variable objects" from 
six OGLE-II LMC "inner" fields (SCI - SC6) covering about 1.3 deg 2 and centered on 
the highest stellar density regions, and ~ 38,000 "variable objects" from the three lowest 
stellar density "outer" fields (SC15, SC19, and SC20). The core of this paper focuses on the 
analysis of the highest density inner fields as a worst case scenario for identifying quasars, 
with the lower stellar density outer LMC fields as a simpler comparison sample. An object 
was defined to be variable if the standard deviation of its light curve from a constant flux 
was at least twice the median photometric error (hereafter called the variability criterion). 
A large fraction of these "variables" are flat light curves with small numbers of outliers, and 



http : //ogledb . astrouw . edu . pl/~ogle/CVS/ 



http : // ogledb . astrouw . edu . pl/^ogle/photdb/ and http : //ogle . astrouw . edu . pi/ ogle2/ dia/ 



http : //ogle . astrouw . edu . pl/ogle2/lmc_ecl/ 
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these are removed by the following cleaning procedure. First, we corre ct the photometric 
errors following the procedure and values from lWyrzykowski et al.l (120091 ). Then we removed 
up to five significant (> 5a) outliers (relative to their adjacent epochs) from each light curve, 
where a typical OGLE-II light curve has 360 epochs. After removing the sources no longer 
satisfying the variability criterion, we are left with 86,301 light curves in the inner LMC 
fields and 25,956 in the outer ones. The inner fields include 584 RR Lyrae (499 "ab" and 
85 "c"), 133 classical Cepheids, 7 double mode Cepheids, 6,745 OSARGs, 374 ellipsoidal 
variable red giants, 46 eclipsing binaries, 699 LSPs, 1,442 LPVs and 15 Be stars. We also 
use the 58 OGLE-II and 984 OGLE-III quasar candidates meeting the variability criterion. 

On further inspection many of these "variable" stars are created by problems in the wings 
of the point spread functions (PSFs) of b right (variable) stars. This has been noted in earlier 
studies (e.g., Ever 2002 : Sumi et "al"1l2005l ). In particular, faint stars near bright variable stars 
tend to vary in phase with the bright source. We apply a three step procedure to remove 
these "ghost" variable objects. First, we mask regions near bright stars. Stars brighter than 
/ »s 12.5 mag are saturated in the OGLE data. These stars have PSF wings extending to 60 
arcsec but they are not formally "detected" in the OGLE catalogs. In order to mask these 
stars, we used stellar objects from the 2 MASS catalog (ICutri et al.ll2003l ) brighter than J = 14 
mag. The masks are defined by rectangular regions with AR.A. = 25, Adecl. = 5 around the 
Two Micron All Sky Survey (2mass) position of each star, where 5 = 120(1 — J/14) arcsec. 
If a variable object is inside one of these regions, it is removed from the list of variables. We 
also masked visually identified regions which showed clear photometric problems (a dense 
clump of variable objects, variables detected on a "bloomed" CCD line, etc.). Our second 
criterion is that no variable could have any other variable object within 6 arcsec. Finally, 
the third criterion is that no variable can have a > 2 mag brighter variable within 15 arcsec. 
These steps could be carried out with greater precision, but cleaning the OGLE catalog of 
spurious variables is not our primary focus even if it is ultimately a limitation and is one 
reason for using our rather conservative variability criterion. One additional problem with 
the OGLE-II photometry is a result of the realuminization of the telescope's mirror. A small 
number of light curves show a step-like change in magnitude between seasons 2 and 3. These 
light curves are identified during the visual inspection we carried out for the fi nal se l ection 
of q uasar candidates. A summary of these photometric problems is also given in lEyerl (120021 ) 
and ISumi et al.l (120051 ) . 
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3. Methodology 



We are interested in modeling random processes leading to variability. Many such pro- 
cesses can be described by the covariance matrix S of the signal, where we will be considering 
the particular case of an exponential covariance 



Si 



V. 



, %3 - o-exp^-ltj - tj\/r) 

between epochs at U and tj used by iKelly et al. J2009h to model quasar light curves base d 
on a standard method fr om the forecasting literature (see, e.g., iBrockwell fc David 120021 ). 
As explained in detail by IKelly et al.l (120091 ). the physics of this process is a random walk 
driving term of amplitude a balanced by a damping time scale r for a return to the mean. 
It is useful to define a 2 = rcr 2 /2 and use r and a as our model parameters - on short time 
scales |At| << r the dispersion between two points is cr| A.t | 1//2 while on long time scales it 
asymptotes to a. Computationally, there is less covariance between r and a than r and a. 
The power spectrum of the process is 



P{f) 



2a 2 r 



2^2 



(2) 



corresponding to white n oise on long time s cales (/ — > 0) and then a oc / 2 fall on short 



time scales (/ > 1/r). IKelly et al.l (120091 ) suggest that r may be related to the ther- 



mal time scale of accret i on dis ks. We estimate thes e parameters using the approach of 



Press. Rybicki Hewittl (Il992al ), its generalizatio n inlRybicki fc Press! (119921 ) and the fast 
computational implementation in lRybicki fc Press! (119941 ) . We will collectively refer to these 
as the PRH method. In the appen dix we summarize this approach, show how it can be used 
to derive the forecasting model of IKelly et al.l (120091 ). and demonstrate that it is more sta- 
tistically powerful than the forecasting approach while still requiring only O(N) operations 
for a light curve with N points. 

Operationally, we fit a light curve by maximizing the likelihood (Equation (IA8I) ) that 
the light curve can be fitted by the process model to determine the parameters r and a (as 
well as their uncertainties, if desired). We include the light curve mean as a simultaneously 
optimized linear parameter (Equations (IA4j) and (1A7I) ). Given the estimates for r and a, 
we have optimal estimates for the mean light curve at both the observed points (Equation 
(1A3P ) and at any other time (Equation (lAllj) ). as well as the variances about these means 
(Equations (1A9j) and (1A12[) . respectively). Where we show these reconstructions and the 
"error snakes" defined by the variances, there are two important points to consider when 
comparing them to the data points. First, these variances are the variances in the mean 
light curve and not the variance of the data relative to the mean light curve. The latter 
quantity, given in Equation (IA10I) . is defined only where there is data and so is ill-suited to 
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showing a continuous light curve. Data points will be scattered relative to the mean light 
curve by the combination of the variance in the mean light curve and the uncertainties in 
the individual data points. Second, the reconstructed light curve is not an example of an 
individual random walk defined by the parameters (r and a), but rather the average of all 
such random walks that are consistent with the observed light curve given its uncertainties. 
The variance in the reconstructed light curve is then the variance of these individual random 
walks about the mean. If we generated individual random walks constrained by the data 
(see PRH), they would track the mean light curve and (statistically) stay inside the "error 
snake" defined by the variances, but they would show more structure on short time scales 
and excursions outside the "error snake" consistent with the estimated variances. 

When we fit the model, we obtain a maximum likelihood estimate of r and a with log 
likelihood hiLt, est . There are two limiting cases for this model. First, as the time scale 
r — > 0, the covariance matrix becomes a diagonal matrix Sij — > cr 2 5ij equivalent to simply 
broadening the photometric errors. We refer to this limit as the "noise" limit with log 
likelihood In L noise . We use a cut of In Lb est > In L noise + 2 to select objects that are modeled 
by the exponential process. Given the variability criterion from Section 2, this cut mainly 
eliminates variable objects with r smaller than the mean epoch separation. The other limit 
t — ► oo, with likelihood lnL^, is the limit where we cannot determine r given the overall 
time span of the survey. This has not proved to be major issue with the OGLE data. We 
cannot distinguish r from r — > oo for only 2% of the ~ 86, 000 OGLE-II variables and 4.5% of 
the OGLE-III quasar candidates. For these objects the time scale r was set to log 10 (r) ~ 5. 

We also check for peri odicity in each ligh t curve using the Lomb-Scargle periodogram 



method as implemented by lPress et al.l (jl992bl ). For each light curve we note the three most 



likely periods Pj (i = 1,2,3) and their probabilities. We used narrow notches at 1, 1/2, 1/3 
and 1/4 of a day to minimize diurnal aliasing problems, while still being able to examine 
RR Lyrae periods. These aliasing solutions could be improved, but this is not the focus 
of this paper. We classify a source as periodic if the estimated period P is shorter than 
200 days and the likelihood for the peak in the periodogram being observed at random is 
log 10 G°p Cr iodic) < ~ 3. We include a limit that periods must be shorter than 200 days because 
many quasars with large r have a "period" satisfying log 10 (p pcri odic) < _ 3 with P > 200 day. 
These solutions corresponds to fitting a sine wave to the data, and as the number of time 
scales t covered by the light curves shrinks, there is an increasing likelihood that a sine wave 
of some period will enormously improve the goodness of fit over no variability, leading to a 
false positive in the periodogram. For periodic variables we should see a strong correlation 
between P and r. 
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4. Quasars 



For our quasar samples we use the 109 known quasars from iKelly et al.l (120091 ) and 
~ 2, 700 mid-IR selected qua s ar can didat es detected in the OGLE-IH LMC and SMC fields 
from iKozlowski &: Kochanekl (120091 ) . The lKellv et al.l (2009) sam ple is a heterogeneous mix- 
ture of MACHO variability selected q uasars ( Geha et ali 2003b and quasars m onitored at 
least in part for reverberation mapping jGiveon et al.lll999l : Peterson et al.ll2004j ). They have 
the advantage of being confirmed quasars, but they are a non-random sample of quasars, 
and many of the reverberation mapping quasars have poorly sampled light curves compared 
to the MACHO or OGLE light curves. The mid-IR selected candid a tes h ave the disad- 
vantage that they are not confirmed quasars. IKozlowski fc Kochanekl (120091 ) classified the 
candidates based on their mid-IR colors (A/B for being inconsistent /consistent with a black 
body), location in the mid-IR color- magnitude diagram (CMD; YSO/QSO for whether the 
object is/is not in the region contaminated with young stellar objects, YSOs) and optical to 
mid-IR color (a/b for whether the optical to mid-IR colors are/are not consistent with spec- 
troscopically selected quasars). Here we only use the "a" sources. Based on the extragalactic 
sources from the AGN and Galaxy Evolution Survey (AGES, C. S. Kochanek et al. 2010, 
in preparation), the QSO-Aa candidates, which are the majority, should be almost entirely 
comprised of quasars, while the other classes are likely dominated by contaminating sources. 

Ideally, we should work with a quasar light curve that is free from contamination by the 
host galaxy. This can be a pr oblem because we fit the light curves in magnitudes rather than 
flux (as did lKelly et al.ll2009l ). In magnitude space, the effects of contamination on the pro- 
cess parameters are non-linear. Extra contaminating flux reduces amplitudes and increases 
time scales, while ove rsubtracting any c o ntamination has the reverse eff ects. For the bright, 
dista nt quasars from iGeha et al.l (120031 ). IKozlowski fc Kochanekl ( 120091 ) and, to a large ex- 
tent, Giveon et al. ( 19991 ). this will be a limited problem. Not only are these very luminous 



quasars, but the ob servations are largely in the r est frame ultraviolet, where host galaxies 



have little flux. The IKozlowski fc Kochanekl (120091 ) sample will also have few sources with sig 



Gorjian et al. 


2008; 


Assef et al. 


2009) 



, MR K 509 JCarone et al1ll996h. NGC 3783 Jstirpe et alill994h . NGC 4051 JPeterson et al. 



20001 ) . and NGC 5548 (jPeterson et aI.ll2002L and references therein), and then either subtract 
the remaining host flux in the spectr oscopic apertures o r add in the host flux from outside 
the spectroscopic apertures based on iBentz et al.l (120091 ). we see significant changes in the 
model parameters. On average, adding the rest of the host flux (subtracting the remain- 
ing host flux in the aperture) reduces (increases) the variability amplitude by an average 
of —0.73 dex (+0.38 dex) and increases (decreases) the characteristic time scale by an av- 
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erage of +0.77 dex (—0.42 dex) for five sources. The effects are larger when including the 
remainder of the host because the spectroscopic apertures used to obtain these light curves 
are already excluding most of the host galaxy flux - the host flux outside the aperture is 
1-10 times larger than the host flux in the aperture. Thus, as with reverberation mapping, 
correlations of the process parameters with other quasar properties will be sensitive to the 
treatment of the host galaxy. We will not consider this problem further since we are focusing 
on the observed properties of quasar variability rather than their detailed interpretation. 

Figure [H shows four exa mples of QSO ligh t curves modeled by the stochastic process. 



These fits are typical, and like lKelly et al.l (120091 ). we find that most quasars are well modeled 



by the process, as illustrated by the x 2 /dof distribution shown in Figure [2j Some care 
is required to interpret this distribution. All objects we consider are likely to be poorly 
described by no variability, so the process model will always produce a reduced x 2 ■ Also keep 
in mind that one limit of the model is simply to broaden the error bars (the "noise" limit), 
which can always allow x 2 /d°f ^ 1 up to the effects of the priors on a and r. Nonetheless, 
the final x 2 /dof distribution is far narrower than the distribution from fitting the linear 
model and few objects are doing so in the "noise limit" corresponding to simply broadening 
the uncertainties. The distribution is somewhat broader than it should be for correctly 
estimated Gaussian uncertainties. For a typical OGLE-III light curve with N = 360 points, 
we expect x 2 /N — 1 ± y2/N = 1 ± 0.07. That the distribution is broader, and relatively 
symmetric, suggests that some of the difference lies in the accuracy of the photometric errors, 
since a 3% misestimation of the uncertainties is enough to cause such a shift. There is some 
skewness to larger x 2 /dof, but this could simply be due our pruning of outliers being too 
conservative - it only takes two 4<r outliers to produce an 0.07 shift in x 2 '/N . Nonetheless, 
we cannot rule out the possibility that the process model is a poor or incomplete model for 
the variability physics in some objects. 

The detection of variability is strongly magnitude-dependent because fainter sources 
require larger amplitudes for variability to be detected. Figure [3] shows the fraction of 
mid-IR quasar candidates that pass the two variability criteria based on either OGLE-II or 
OGLE-III data. One advantage of the process model is that we can explicitly calculate the 
completeness of a survey as a function of the parameters using Monte Carlo simulations. At 
each point on a grid of the process parameters (a and r), we randomly generate N light 
curves using the temporal sampling of a randomly selected real light curve. Generating 
simulated light curves of zero mean is trivial, since the next point 

s i+1 = Sl e~^ + G [a 2 (1 - e- 2A ^)] , (3) 

where At is the time interval and G(x 2 ) is a Gaussian deviate of dispersion x. The observed 
light curve is then Hi = Si + G(nf) where rii is the observational noise. The chain is initiated 
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by si = G(a 2 ). We use the mean measurement error at a given magnitude, but one could use 
the error estimates from a light curve matched in magnitude to the trial. Next, we apply our 
variability criterion, that the rms of the light curve must be twice the measurement errors, 
and that \iaL best > \B.L noise + 2, to each trial light curve. The completeness is simply the 
fraction of the trials which pass both selection criteria. For speed, we estimate the second 
criterion by using the likelihood at the input parameters for Lb est and the likelihood for the 
model with the same asymptotic variance but r — > for L noise . 

Figure H] shows the resulting completeness estimates for OGLE-II and OGLE-III light 
curves at / = 16.5, 17.5, 18.5, and, 19.5 mag, where the typical photometric uncertainties 
(including the rescalings of Section 2) are a p hot = 0.010, 0.022, 0.050, and 0.12 mag for OGLE- 
II and a phot = 0.008, 0.017, 0.037, and 0.087 mag for OGLE-III. The general structure of 
the completeness limit is easily understood. The variability criterion can be approximated 
by noting that in the continuum limit, where we integrate over a continuous light curve 
rather examining discreet samples, the variance of a light curve described by an exponential 
covariance matrix relative to its mean is 



Var (c(t)) = a 



2 2 

1 — - + 3 (1 - exp(-s)) 



x ar 



(4) 



where x = At/r is the ratio of the survey duration to the time scale r. When the survey 
duration is long compared to the time scale (x ^> 1), the signal variance is simply cr 2 , 
and light curves above the diagonal line with a 2 + o~ 2 phot > ^o~ 2 hot will satisfy the variability 
criterion. When the survey is short (x <C 1), the observed variance is limited by the survey 
duration to <t 2 x/3 = cr 2 At/6 and the completeness is high if a is sufficiently large. The 
criterion that the process is distinguishable from simply increasing the error bars eliminates 
sources in the corner with short time scales and high amplitudes. Note that the variability 
criterion is very conservative, because its selection limits do not improve with additional 
data. Using only the second criterion, based on whether the process model fits better than 
simply expanding the measurement errors, would lead to higher completeness for longer r 
and smaller a but at the expense of higher false positive rates due to systematic errors in the 
data. Compared to the OGLE-II survey, the OGLE-III survey is more sensitive and covers a 
longer temporal baseline but with a larger epoch spacing, so it is more complete for large r 
and smaller a, while less complete for small r. The completeness limit is largely determined 
by the variability criterion except for the corner with high amplitudes and short time scales. 
While not shown, the distribution of variable sources as a function of magnitude follows 
these completeness limits closely. 

Figure [5] shows the distribution of quasars in an optical CMD, 07 — r parameter space 
(shown in detail in Figure [6]), 07 — magnitude space and, where available, the relationship 
between 07 and ay at fixed r. Because the V-band light curves are poorly sampled compared 
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to the /-band, we determine <5y using the r derived for the /-band light curve. To better 
understand the contamination problems we discuss later, we s uperpose a CMP of LMC stars 
based on the HST Local Group Stellar Photometry archiv^ ( Holtzman. Afonso &: Dolphin 



20061 ). On the CMD one can quickly identify the red clump (RC), the red giant branch (RGB) 



and the top of the main sequence. Since these are derived from small Hubble Space Telescope 
(HST) fields, the sequences of higher luminosity stars are not significantly populated. 

We see that quasars occupy a well defined locus in the a-r space, lying in a band with 
— 21og 10 (<j) + 0.3 < log 10 (r) < — 21og 10 (cr) + 1.6 that roughly corresponds to an asymptotic 
variability of a ~ 0.1 mag. Variability time scales longer than a few 1000 days are rare, 
although r is indeterminate but long for about 5% of the sample. The estimate of r becomes 
very uncertain as it approaches the survey length. At very short r, systems are lost because 
they fail the lnLj, est > lnL noise + 2 criterion. If we focus on the QSO-A sources, which are 
overwhelmingly going to be real quasars, we can define a shaded region in Figure [5] (Cut 2 
in Section [6]) that encompasses most of these sources (69% of those identified as variable in 
OGLE-III). The more contaminated QSO-B and YSO candidates are less likely to fall into 
this region (see Figure El Tables [H and [2]) , providing added evidence of their higher levels of 
contamination by LMC sources. The mid-IR quasar candidates tend to lie in the Hertzsprung 
Gap of the CMD, consistent with the typical color 0.4 < (V — I) < 1.0 mag of z < 3 quasars. 
They also show relatively achromatic variability between the V and / bands, with more 
chromaticity for small a and larger r sources. Tables [T] and [2] summarize the properties of 
the mid-IR selected quasars for both OGLE-II and OGLE-III light curves. If we compare 
the locations of the quasars in the space of the variability parameters (Figure ED to the 
completeness estimation of Figure HJ it is easy to understand the magnitude dependence of 
the fraction of variable quasars in Figure [31 

The distribution of the quasars is a combination of intrinsic properties (variability at 
rest wavelength A for a black hole of some mass and luminosity), distance (magnitude and 
scalings between rest-frame and observed-frame properties), and selection effects (Figure H]). 
The observed magnitude of the quasar depends on distance and /("-corrections, and the 
variability parame ters have redshift corrections (r rest = r/(l + z) and a rest = a(l + z) l l 2 
(jKelly et al.ll2009l )) and potentially ZT-corrections for any scalin gs of the parameters wit h 
wavelength, such as those seen in average structure functions (e.g. JVanden Berk et al.ll2004j ). 



5 http://ganymede . nmsu.edu/lgj - we selected two HS T fields lmc_u2xq05 and lmc_u2xq06, each 
observed with the WFPC2 for 1540 and 1860 seconds in F555V and F814W, respectively. Since in both 
cases the PCI chip was pointed at a globular cluster, we used the data from the three remaining WF chips, 
transformed to V- and /-band. These data were then binned into 0.05 mag bins in both magnitude and 
color and are shown as smoothed contours of 1, 5, 10 and 20 objects per bin in Figures O [7] and [5] 
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Our incompleteness depends on observed rather than intrinsic properties, and this makes it 
highly likely that our estimate of the region occupied by quasars is largely correct despite 
the increasing problems with completeness at fainter magnitudes. In essence, the majority 
of quasars which we lose due to incompleteness correspond to quasars with similar intrinsic 
properties at lower redshifts where they are bright enough to be included, rather than quasars 
with intrinsically different properties. The numbers of quasars observed to have a given set 
of parameters are, however, strongly dependent on the selection function, and reconstructing 
the true distribution will require redshifts and (essentially) a V/V max method applied to the 
variability selection procedure. A significant advantage of our approach is that such problems 
can be quantitatively addressed using the variability model. The LMC sources we consider 
in the next section are an illustration of the opposite limit from the quasars. Since all LMC 
sources are at the same distance, we either detect them in this variability space given their 
magnitude or not. 



We can also examine the distribution of the iKelly et al.l (120091 ) quasars relative to the 
mid-infrared sample in Figure [6j They generally lie in the same region (79% are in the 
selection box), but their distribution in that region is very different. Most of the MACHO and 
PG quasars lie at the large r end of the distributio n, with a few of the reverberation mapping 



quasars dribbling toward short r. The MACHO (jGeha et al.l 120031 ) quasars are variability 
selected, but the selection criterion may be a biased one compared to a (roughly) magnitude 
limited sampling of quasars. Figure [6] also shows the effects of adding or subtracting host 
galaxy discussed above. 



5. General Variability of Sources in the LMC 

In Section HI we demonstrated that quasars occupy a limited region of the variability 
parameter space. In order to select quasars using variability, we need to examine how other 
variable sources populate these parameter spaces. In Figure we show the location of 
the general variable populations from the inner LMC fields for the same combinations of 
variables used for quasars in Figure [51 The distributions of outer field variables are similar 
but of lower density. Where possible we have identified the variable type using the catalogs 
discussed in Section [2j Tables [1] and [2] summarize the statistics of selecting sources in 
these fields. It is critical to clean the input catalogs, since the variability criterion identifies 
~ 86, 000 inner field variables before we apply the bright star masks, defect masks, and the 
isolation criteria described in Section (2J while there are only ~ 13, 000 left afterward. Of 
these, ~ 65, 000 (no masks) and ~ 10, 000 (with masks) pass the restriction that the process 
model fits the data significantly better than simply expanding the photometric uncertainties 
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(In Lbest > \n-L noise + 2). While the tables show the statistics for both cases, we only discuss 
the results for the sources including the masking and isolation criteria. 

In Figure we can quickly identify a number of sequences or overdensities of variable 
objects. We separate them, for clarity, in the four panels of Figures I8HT21 In these figures 
we plot rare samples last so that they are not masked by the common ones. Quasars are 
shown in multiple panels to maximize their visibility compared to other populations, and 
Figures and are always available for reference. In the top-left panel of each figure we 
show two samples of OSARGs, where "a" corresponds to th e asymptotic giant branch (AGB) 



objects and "b" to the RGB objects (jSoszyriski et al.l 120071 ). They are approximately 3 mag 



bright er than the RC and also r edder. In this panel, we also show slowly pulsating B stars 



fSPB. iKolaczkowski et al.ll2006l ). active giants, active subgiants, and RS CVn binaries (AG 
and ASG, Drake 20061 ). There is no formal OGLE catalog of these variables, so we defined 



them based on the selection boxes shown in Figure [H The OSARGs share the same area 
of the CMD with the LSPs, ELLs (top-right panel), and LPVs (bottom-right panel). The 
OGLE eclipsing binary samples are located just above the RC, and, slightly higher, one finds 
ellipsoidal variable red giants (top right). Along the instability strip with an approximate 
color of (V — I) = 0.5, one can find Cepheids, ~ 3 mag brighter than RC, and RR Lyrae, 
which are slightly fainter then RCGs. The fundamental-mode Cepheids (RR Lyrae), denoted 
as "FU" f a b"), and first-overtone ones, denoted as "FO" ("c"), form distinct sequences. 



Geha et al.l (120031 ) in their attempt to select QSOs based on variability found that Be/Ae 
stars are the main source of contamination. Be stars can be found in the bottom-right panel 
just above the main sequence. In the bottom-left panel we present all the remaining variable 
objects. These should be quasars, other unclassified variables and sources with unidentified 
systematic problems. For comparison, we also show the OGLE-III quasar candidates in three 
of the panels, excluding the panel of unclassified sources. 

The top-right panels of Figures [3] and [7| show the positions of objects in the "modified 
amplitude" -magnitude (a— I) space. Once again, a number of source groupings are easily 
noticed, where the most prominent ones are OSARGs, LSPs, LPVs, RR Lyrae, and Cepheids 
(color coded in Figure ITOl . OSARGs, ELLs, LSPs, and LPVs overlap, although they can be 
partially distinguished based on their periods (Figure |T2j). In a— I space Cepheids and RR 
Lyrae occupy distinct areas. 

We show the positions of all the variable stars and quasars in the space of the stochastic 
model parameters in Figures [5] and [7| (bottom- left panel). In this space of modified "char- 
acteristic time" and "variability amplitude" (r — a), we color coded all variability classes 
and present them in Figure [91 A majority of OSARGs and a substantial fraction of ELLs, 
LSPs, and LPVs can be easily distinguished from the quasar locus, although there is heavy 
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contamination for long r. Note also that this is also the region occupied by most MACHO 
quasars, and this may be one cause of low yield of variability selected quasars in Geha et 
al. (2003). Cepheids, RR Lyrae, some LPVs, Be stars, ECLs, SPBs, and (3 Cep, AGs and 
ASGs also overlap the quasar locus area. Note that fundamental mode and first overtone 
Cepheids and RR Lyrae can be distinguished from each other in a — r space. 

The bottom-right panels of Figures [5] and [7| show the variability ratio between the V- 
and /-bands as a function of /-band variability. The majority of LMC objects show higher 
variability amplitudes in the V-band than in the /-band, as we would expect since stellar 
variability is dominated by temperature changes. Quasar variability is less chromatic, so the 
amplitude ratio is another tool for separating quasars from variable stars. For these LMC 
fields, this is not critical because of the large magnitude differences between many of the 
variable stars and quasars. However, in a field dominated by Galactic sources spread over a 
wide range of distances, amplitude ratios have the advantage of being distance independent. 
The color-coded variability classes are shown in Figure [TT1 

Figure [T2] shows all the variable sources and quasars (blue dots) in the space of the 
characteristic time-scale r and the most likely period P derived from a standard periodogram. 
As noted in Section [3X quasars never have significant probabilities of being periodic sources 
unless the "period" is long, P > 200 days. On inspection, these light curves have fluctuations 
on long time scales that are quasi-periodic over the extent of the light curve (particularly 
with yearly or half-yearly aliases), so a sine wave fit to the data is a significantly better 
fit to the data than a constant. Even the quasars with falsely probable periods show no 
correlation between the period and r. This leads us to define a periodic source as one with 
P < 200 days and log 10 (p pe riodic) < — 3. On short time scales (r or P < 1 day) we see no 
real correlations because the periodograms have diurnal aliasing problems and many sources 
fail to pass the In Lf, est > In L noise + 2 variability selection criterion. For periods of days to 
tens of days, many sources follow a scaling with r oc P 2 , which is surprising since we expect 
t oc P by dimensional analysis. 

Our one hypothesis for this scaling is that it can be produced if the exponential process 
model is dominated by how well it fits the peak of the correlation function. Near the peak, the 
auto-correlation function of a sine wave is proportional to At 2 /P 2 + const, while that of the 
exponential covariance matrix is At/r + const, potentially driving the r oc P 2 scaling. There 
is a regime where the longer period Cepheids have r oc P, and some of the multiperiodic 
sources (OSARGs) appear in multiple groupings. Monte Carlo experiments fitting sine wave 
light curves show that the relation between r and P is strongly affected by the temporal 
sampling, with the r oc P 2 region shifting to longer periods as the mean epoch spacing 
increases. Despite some effort we could not derive an analytical model for the scaling. It 
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is clear, however, that the relationship between r and P depends on sampling rather than 
being universal. 



6. Variability Selection of Quasars 

We can now combine the results from Sections H] and to define selection criteria for 
isolating quasars from other variable sources. Presently this is limited to searching for candi- 
date quasars in the 1.3 deg 2 (0.65 deg 2 ) covered by the OGLE-II inner (outer) fields in order 
to explore how well we can find quasars while avoiding contamination. We limit the analysis 
to sources with / < 19.5 based on the completeness estimates in Section |H There should be 
very few quasars in this region, so this is really a test of contamination levels in a field with 



very high densities of stellar variables. Based on iRichards et al.l (120061 ) we expect approxi 



mately 1, 6, and 20 quasars per square degree with / < 17.5, 18.5, an d 19.5 mag respectively. 



which agrees well with the statistics of the mid-IR selected quasars (IKozlowski fc Kochanek 



20091 ). We will simultaneously apply the selection method to the broader sample of OGLE- 
II and OGLE-III quasar light curves to examine their effects on completeness. We select 
quasars in four steps, starting with two that will be generic to any set of light curves from 
any location (Cut 1 - non-periodicity, Cut 2 - variability properties of quasars), and then 
applying those that depend on the specific field (Cut 3) and the availability of light curves 
in multiple bands (Cut 4). We will not apply a color selection criterion, since one goal of 
any wide application of this method would be to find q uasars at redshifts w here their colors 



are crossing the stellar color distribution near z ~ 2.6 (IRichards et al.l 120061 ) . 



Cut 1. We start by eliminating periodic sources. A periodic source is defined to have a 
maximum likelihood period P < 200 days with ln(p periodic ) > —3. This removes 76% of the 
sources, including, 99% of RR Lyrae, 100% of Cepheids, 91% of ELLs, and 100% of ECLs. 
It does not remove 20% of OSARGs, 69% of LSPs, and 36% of LPVs. 

Cut 2. Next we isolate sources with the variability properties of quasars, defined by 
the region in a — r space bounded by — 21og 10 (<3") + 0.3 < log 10 (r) < — 21og 10 (cr) + 1.6, 
log 10 (<r) > —1.1 and r > 2 days. This quasar region is shown in gray in Figures IU O [61 
[7] and [9j The lower limits are designed to remove most Cepheids, RR Lyrae, OSARGs, 
ellipsoidal variable red giants, long secondary period variables, slowly pulsating B stars, (3 
Cephei stars, active giants and subgiants, and YSO-(AB) objects. The upper limits remove 
many of the long period and long secondary period variables. This criterion will still leave 
significant contamination from several of the variable star populations along the left (low 
a) and bottom edges. These are blue variables, active giants and subgiants and RR Lyrae. 
This cut removes 97% of ELLs, 92% of OSARGs, 45% of LSPs, and 41% of LPVs passing 
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Cut 1. 

Cut 3. The third cut takes advantage of the fact that almost all the contaminating 
sources are at a common distance and largely consist of relatively bright giant stars. It is 
defined by a wedge in <r— magnitude space that isolates the OGLE-III quasars. We first 
require quasar candidates to be fainter than / = 16 mag since such bright quasars are very 
rare, while the OSARGs, classical Cepheids, eclipsing binaries, ellipsoidal variable red giants 
and long period variables are all brighter than / = 16 mag given the distance to the LMC. A 
second limit, I > 3\og w (a) + 18.5, separates the RR Lyrae, double mode Cepheids, ECLs, 
some ELLs, SPBs, (3 Cep, AGs and ASGs, which are to the right of the quasar region. As 
in Cut 1, we use the same lower limit on the scaled amplitude of log 10 (<r) > —1.1, which 
partially eliminates Be stars. This cut removes 100% of the ELLs, LPVs, LSPs and OSARGs 
passing Cut 2. Cut 3 would need to be adjusted for other background environments, such 
as the Galactic bulge, where similar variable stars are ~ 4 mag brighter than in the LMC 
(depending on the amount of extinction). 

Cut 4- For the objects with well determined F-band amplitudes <5y we can also separate 
quasars using the ratio of the V- and /-band amplitudes. We define an area in (5y/<7j)-07 
space with —1.1 < log 10 (oj) < 0.6 and log 10 (5y/aj) < — 0.2 log 10 (<3j) + 0.2, where we select 
objects with a smaller ratio of the V- to /-band variability amplitudes than is typical for stars 
(Figure [TU gray area). This cut removes most Cepheids and RR Lyrae and a substantial 
fraction of other common variable stars. 

The effects of the selection functions on both the LMC variable population and these 
quasar samples are summarized in Tables [1] and [2] both for their effects when applied se- 
quentially (Table [T]) and individually (Table [2]). We also summarize the effects on the sub- 
categories of the mid-IR selected samples, where we expect to eliminate more of the heavily 
contaminated QSO-B and YSO classes than the relatively high purity QSO-A class. Recall 
from Section H] that we lose a large fraction of the quasars from the two requirements for 
being considered a variable because most quasars are relatively faint compared to the depths 
of the OGLE data. 

The first test of a selection method is its completeness. We evaluate this based on the on 



the ~ 2, 700 (~ 300) mid-IR selected quasar candidates from iKozlowski fc Kochanekl (120091 ) 
with OGLE-III (OGLE-II) light curves as well as the smaller iKellv et al.l J2009h sample. 
The biggest problem for variability selected quasars from these samples is that most quasars 
will be faint. While we can identify 2716 (310) of our mid-IR quasars in the OGLE-III 
(OGLE-II) data, only 1139 (215) are bright enough (/ < 19.5 mag) to plausibly detect 
quasar variability (see Figure [3]). For the full OGLE-III (OGLE-II) samples, 984 (58) pass 
the basic variability criterion of Section [21 while 635 (58) of the / < 19.5 mag candidates 
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pass this criterion. Adding the noise criterion, lnL^t > lnL no j se + 2, we find 619 (848) 
for OGLE-III and 55 (55) for OGLE-II, brighter than / < 19.5 mag (no magnitude limit). 
The difference in the fraction of variable objects in the two phases of the OGLE survey is 
explained by the added depth (0.3 mag) and duration (4 versus 8 years) of the OGLE-III 
survey. By design, the remaining cuts have little effect on the quasar candidates, although 
we lose a larger fraction of the non-QSO-A sub-samples, as we would expect given their 
higher levels of contamination. After applying the first three cuts to the variable objects, we 



are lef t with 65% of the QSO-A sample. We can apply only Cuts 1 and 2 to the iKelly et al. 



(120091 ) sample, and 65% of the sources pass the cuts. Cut 4 removes an additional 7% of the 
OGLE-III QSO-Aa objects. Finally, all four cuts are met by 61% (63%) of the OGLE-III 
QSO-Aa quasar candidates with no faint magnitude limit (within 16 < / < 19.5 mag) and 
17% of the OGLE-II candidates. Thus, following this approach can yield variability-selected 
samples of quasars with moderate completeness. 

Next we examine how well we reject other variable sources, focusing on our analysis of 
the central region of the LMC. We only discuss the sample including masks and isolation 
criteria (see Section [2]), but report the results for the full, unpruned sample as well in Tables [T] 
and [2J The selection cuts are very effective at reducing the numbers of candidate sources, 
going from ~ 10, 000 (~ 600) sources passing the variability criterion in the inner (outer) 
fields to only 58 (2) candidate quasars. In the outer fields with fewer stars, the contamination 
is remarkably low. 

Of the 58 remaining OGLE-II candidates in the inner fields, 3 are also mid-IR selected 
quasars. There are 63 (41 QSO-A) OGLE-II sources in the same fields, 24 (12) passing the 
two initial variability criteria, showing that the biggest problem for quasars is recognizing 
variability given the depth of the data. The mid-IR quasar selection method is not perfect, 
in particular it tends to miss lower redshift quasars with significant mid-IR emission from the 
host galaxy ( Gorjian et al. 2008 ; Assef et al. 2009^. so some of the remaining, unidentified 



candidates may be real quasars. For the iRichards et all (120061 ) quasar number counts, we 



expect of order 1, 8, and 25 quasars brighter than / < 17.5, 18.5 and, 19.5 mag in this 
1.3 deg 2 region, compared to 6, 16 and 41 QSO-Aa candidates in OGLE-II. We examined 
the light curves of all 58 objects that passed the four criteria. On visual inspection, the 
majority seem to be created by photometric problems in the wings of other stars, long 
period variable stars ju st above the RCG, stars in the wings of high proper motion stars (see 



Ever fc Wozniakl 120011 ). blue and Be stars, eruptive variables, novae /super novae, nova-like 
variable, unidentified types of variable stars, and quasar candidates. Of the non-mid-IR 
sources, only one seems likely to be a quasar, the OGLE-II (OGLE-III) object LMC_SC5 
253536 (LMC167.2 36755, at J2000 05:24:24.17 -69:55:56.9). It has a brightness of / = 18.53 
(18.37) mag in OGLE-II (OGLE-III). Its color of (V - I) = 1.38 (1.06) mag is somewhat 
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redder than expected for z < 3 quasars. We show both the OGLE-II and OGLE-III light 
curves in Figure [131 

In the outer OGLE-II fields we are left with only 2 candidates. On inspection we find 
no plausible QSO candidates. None of the OGLE-II light curves of the 43 mid-IR selected 
candidates (29 QSO-Aa) from the outer fields passes all the criteria. In fact, only two 
QSO-Aa and two YSO-Aa objects with OGLE-II light curves pass the two initial variability 
criteria. 



7. Summary 



Periodograms to recognize periodically varying sources and begin their classification are 
a standard astronomical tool. Here we introduce an equally simple approach to classifying 
continuously varying sources that are not periodic and explored the properties of large sam- 
ples of variable sources in the LMC and SMC using data from the OGLE-II and OGLE-III 
microlensing projects. We modeled the sources as a damped random walk, a stochastic pro- 
cess with an exponential covariance matrix characterized by a scaled amplitude a, a time 
scale r, and an asymptotic var iance on long time scales of a = a{r/2) 1 / 2 . This model was in- 
troduced by lKelly et al.l (120091 ). who used a forecasting approach to estimate the parameters 
that is not as statistically powerful as the PRH method used here. Extracting the damped 
random walk parameters is no more complex than obtaining a periodogram, so it can be 
easily used to characterize any light curve. 



We first applied this method to a much l arger sample of quasars thanlKellv et al.l (120091 ) 
by using the mid-IR selected quasars from Kozlowski fc Kochanek ( 20091) that also have 
OGLE-II or OGLE-III light curves. We have verified the suggestion by iKelly et al.l (120091 ) 
that quasar light curves are well modeled by the damped random walk process. Moreover, we 
find that quasars occupy a well defined region of the parameter space of the model, allowing us 
to design a simple quantitative method for the variability selection of quasars. The details of 
the approach will likely require some modifications as larger samples of quasars are analyzed 
and the method is applied in regions with different sources of stellar contamination, but the 
broad outline should survive. We also note that the model allows quantitative estimation 
of completeness in variability selected quasar samples and can estimate the detection ranges 
needed for computing luminosity functions or other population statistics. 



While IKelly et al.l (120091 ) are certainly correct that a damped random walk correspond- 
ing to an exponential covariance matrix models quasar variability well, it is less clear that 
the correlations between the model parameters and the luminosity or black hole masses of 
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quasars are robust. Since we lack redshifts for our quasar sa mple we cannot ex plore this 
in detail, but we note two p otential systematic problems in the Kelly et all (120091 ) analysis. 
First, the distribution of the iKelly et all (120091 ) quasar sample may be a biased sampling of 
the quasar distribution. It appears to be biased toward sources with long time variability 
time scales compared to the mid-IR quasar sample. Second, most of the leverage for deter- 
mining the slopes in the relations between r and a and the quasar mass and luminosity in 
Kelly et all (120091 ) come from the relatively low black hole mass and luminosity r everberation 



mapp ing quasars. Host galaxy contamination is a major issue for these systems ( iBentz et al. 



20091 ). and we find that this also holds for estimates of the variability process parameters 
because light curves in magnit udes are non-linearly distorted by contamination. Thus, the 
slopes of the IKelly et all (120091 ) correlations may strongly depend on the treatment of these 
hosts. Just as for the reverberation mapping studies, we should really measure the variabil- 
ity parameters of the quasar with no diluting contribution from the host galaxy, but this is 
virtually impossible to do for large samples. This problem could be avoided by modeling 
the time variable flux rather than modeling the time varying magnitudes, but the meaning 
of a then becomes strongly distance dependent and fractional variability, as represented by 
variability in magnitudes, probably has greater physical meaning. These issues are not a 
concern for simply classifying light curves. 

We developed a set of selection cuts to try to identify quasars based on their variability, 
and 63% (17%) of the variable quasars in the QSO-A sources in the OGLE-III (OGLE-II) 
surveys pass the selection criteria. In a lower stellar density field this can be improved 
by using more liberal selection criteria. If we apply the same cuts to the 10,000 (600) 
OGLE-II variable sources in these high density inner (lower density outer) fields, only 58 
(2) candidates survive, of which 3 (0) are mid-IR selected quasar candidates. After visual 
inspection, only 1 (0) of the other sources is a quasar candidate. Most of the rest are artifacts 
that could be avoided by better pruning the input catalogs. The remaining sources are 
blue and eruptive variables, variable giants and supergiants. While 85-93% contamination 
(depending on how we count remaining artifacts) seems poor, it was achieved when searching 
for rare extragalactic sources in an extremely high density stellar field (~ 10 6 stars and ~ 
8,000 variable stars/deg 2 compared to ~ 20 quasars/deg 2 with 16 < / < 19.5 mag), and 
the bulk of the false positives were easy to recognize by visual inspection. If we add a color 
cut for z < 3 quasars (0.4 < (V — I) < 1.0 mag) we can also greatly reduce the residual 
contamination. In the outer LMC fields (~ 10 5 stars and ~ 10 3 variable stars/deg 2 ) or an 
extragalactic field (> 10 3 stars and a small number of variable stars/deg 2 ) there will be little 
difficulty with contamination. Moreover, it is clear that the OGLE-II data are too shallow 
to search for quasars in small areas compared to OGLE-III. 



Clearly the next steps are to apply this approach systematically. Here we carried out 
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a limited analysis of OGLE-II and OGLE-III data pending the general availability of the 
full sample of OGLE-III light curves. We are also in the process of obtaining redshifts for 
many of the mid-IR quasar candidates. An obvious next step is to characterize the entire 
OGLE-III survey. In C. L. MacLeod et al. (2010, in preparation), we analyze the nearly 
9,000 known quas ars with light curves from the SDSS to explore the correlations found by 



Kelly et al.l (120091 ) . These SDSS data can also be used to explore the problems of background 



sources and their classification for typical extragalactic fields rather than the center of the 
LMC. 
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A. The PRH Method 



The basic idea flPress. Rvbicki fc Hewittlll992al : iRvbicki fc Presslll992l . 1 19941 ) is to sup- 
pose that we have data y that is due to an underlying, true, signal s, measurement noise 
at each point n, and a general trend defined by a response matrix L and a set of linear 
coefficients q, thus, y = s + n + Lq. For example, we will use the linear coefficients to 
optimally remove the light curve mean, so we have one linear coefficient q\ for the mean, and 
the response matrix is simply a column vector Ln = 1 with an entry for each data point, 
i = 1-..JV0 The intrinsic variability has covariance matrix S = (ss) and the noise has 
covariance matrix N = (nn). We will be explicitly relying on Gaussian statistics in order 
to add normalizing pre-factors to the squared difference statistics focused on by PRH. By 
definition, we know that 



P(s) oc \S\~ 1/2 exp 



s T S~h 



and that -P( n ) oc \N\ 1 exp 



n 



(Al) 



6 Additional uses of the linear parameters follow trivially. To remove a linear trend, q\ +q2(t— to), from the 
light curve, we simply add a second column = U — to to the response matrix. If we have two light curves 
with a possible offset, we could use separate means q\ and q 2 for each segment, and then (Lii,L i2 ) — (1, 0) 
for data in segment 1 and [L%i, La) = (0, 1) for data in segment 2. 
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Thus, the probability of the data given the linear coefficients q, the intrinsic light curve s, 
and any other parameters of the model p (e.g. r and a) is 

/ r / >s, T ^- l >i 4- n T A/" -1 n\ 
P (y|q, s, p) oc ISiVr 1 / 2 y d»n6 (y - (s + n + Lq)) exp ^ 1 J . (A2) 

After evaluating the Dirac delta function, we complete the squares in the exponential with 
respect to both the unknown intrinsic source variability s and the linear coefficients q. This 
exercise determines our best estimate for the mean light curve, 

s = 1 SC- 1 (y-Lq), (A3) 

and the linear coefficients, 

q = (L T C- 1 L)- 1 L T C- 1 y = C^C^y, (A4) 

where C = S + iV is the overall covariance matrix of the data and C q = (L T C _1 L) _1 . With 
these definitions we can factor the argument of the exponential into 

P(y|q,s,p) oc \SN\ exp I j- ^ — I (A5) 

where 

Cl l = C- 1 - C^LCg^C- 1 (A6) 

is the component of C that is orthogonal to the fitted linear functions, the variances in the 
linear parameters are 

(Aq 2 ) = (L^L)- 1 = C v (A7) 

As = s — s and Aq = q — q. We can marginalize the probability over the light curve s and 
the linear parameters q under the assumption of uniform priors for these variables to find 
that 

P(y|p) oc \S + N\-^ |L^L|- 1/2 exp (-^^) , (A8) 

where for the exponential model the remaining parameters p are r and a. This is the 
likelihood we optimize to determine r and a. We also know the variances in the estimate 
for the mean light curve 

(As 2 ) = S - S T C ± S (A9) 
and the variance between the data and the estimated light curve is 

<(y - s - Lq) 2 ) = N(l - C^C^C^il - C^C'^N. (A10) 
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Our only addition here compared to PRH is keeping track of the normalizing prefactor of 
the exponential. 



The relation to [Kelly et al.l (120091 ) comes from using this formalism to predict the value 
of the time series at an unmeasured time. The simplest means of doing so is simply to pad the 
data vector with additional fake points y/ that have infinite measurement uncertainties 
in the sense that iV -1 — > for these points. For simplicity we consider the case without 
additional linear parameters. We partition the signal and noise matrices as S dd , Sff, S d f and 
Sdf and , N7f = 0, N^ 1 = and N^ 1 = for the data-data, fake-fake, data-fake, and 



fake-data blocks of the matrices. Substituting into Equation (1A3I) . we find that the estimate 
of the true light curve for the measured data points is the result we would obtain without 
having padded the data vector, and that the estimate of the light curve at the additional 
points is 

s f = S fd (S dd + N dd y l y d (All) 
with variance relative to the mean light curve of 

(As}) = S ff - S fd (S dd + N M y l S df . (A12) 

These two expressions define the mean light curves and the "error snakes" shown in Figures [T] 
and [T3J Mathematically, the mean light curve is the weighted average of all process light 
curves described by parameters p that are statistically consistent with the data, and the 
variance is the scatter of these light curves about this mean. 

We can use these results to "forecast" the light curve by noting that we can express the 
expected value at epoch % + 1 in terms of the measurement at i and the forecast for i based 
on the i — 1 earlier data points. In component language, 



where j, k = 1 • • • i — 1, and 



Sij(S + N)j£vk (A13) 



s i+1 = S i+lj (S + N)j k l y k (A14) 



where j, k = For an exponential covariance matrix, Si+ij = aSij where a = 

exp(— \U + i —ti\/r). If you partition the S + N matrix of Equation flA14j) into the i — 1 xi — 1 



part used to estimate Sj and the remaining column vectors and scalar and then use the 
standard method for inverting a partitioned matrix you find that 

a(Su -Ti) 

s i+1 = asi + —{Vi - Si), (A15) 

where Tj = S T (S + N)~ l S an d Oj is the measurement error for y^. To complete the identifi- 



cation with lKelly et al.l (120091 ). we note that the variance between the data point yi and the 
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prediction Sj is 



((Vi 



2 ) = o 2 + -Ti = a 2 + ^ 



(A16) 



so 



Kellv et all (120091 ) 's = C« - IV Note that Equations flATUD and flAT6|) differ because 



Sj in Equation (1A13I) does not depend on the data point yi, but it does depend on that data 
point in Equatio n flAlOj) . Making this substitution, we recover the forecasting expression 
from iKelly et al.l (120091 ) where 



s i+1 = asi + „ - Si 

af + 

Similarly one can work out the variance in s i+1 to find that 



(A17) 



t+i 



a 



ct, 2 + 



(A18) 



which then gives the second of the forecasting equations in IKelly et al.l (120091 ) 



n 



i+l 



a 



+ a 2 ^ 1 



(A19) 



where Qq = Cu is the diagonal element of the covariance matrix. IKelly et al.l (120091 ) initial- 
ized their predicted time series at the process mean, while the PRH approach would initialize 
it to the first measurement. For this forecasting approach, the likelihood of the data given 
any parameters is then 



P(y|p) ocII^ + ^-^exp 



(Vi ~ Sj) 2 
2(a t + a 2 ) 



(A20) 



instead of Equation (1A8I) . 



Having demonstrated that the IKelly et al.l (120091 ) forecasting method can be derived 
from the PRH approach, we next address the relative virtues of the two methods. An 
apparent advantage of the forecasting approach is that it is computationally O(Ndata) to 
compute Equation (IA20I) . while at first glance it is 0(N^ ata ) to compute Equation (1A8R be- 
cause of the matrix inversion. This is not the ca se for the exponential correlation function 
(Equation ([T])) because its inverse is tridiagonal (jRybicki fc Press! I1994J ). If we rewrite the 
matrix (S + A^) _1 in Equation (1A20j) as S~ 1 (S~ 1 + A^ 1 )^ -1 , then for a diagonal (or even 
tridiagonal) noise matrix N, all the matrix calculations, including the determinant, in Equa- 
tion (1A8jl . can be carried out in O(Ndata) operations. The total number of operations and 
the complexity of the implementation is somewhat higher if one includes the steps needed 
to automatically marginalize over the light curve mean using the linear parameters, but this 
then avoids having it as a parameter which must be optimized as part of the fit. Using 
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the fast method from iRybicki fc Press! (Il994j ) requires the data to be time ordered and non- 
overlapping (ti < i 2 < t 3 • • • < tjy) since S is singular if any £j > tj for i < j - data points at 
identical (or nearly identical, in the sense that \ti — tj\/r << 1) times should be combined 
before applying the analysis, and some care is required to stably compute the determinants. 

Balancing the computational complexity is that the PRH method is more statistically 
powerful because it uses all the information in the light curve simultaneously. In essence, for 
parameter estimation we should use all available data to predict s i+ i not just the preceding 
data. We can demonstrate this by generating Monte Carlo light curves using the exponential 
process, estimating the parameters using both methods and examining the relative scatters 
in the recovered values. We use priors of P(r) — 1/r and P(a) = 1/cr, as logarithmic priors 
are generally standard for positive definite but otherwise scale free variables. This differs 



from lKelly et al.l (120091 ). who assumed a uniform prior for a and a = exp(— At/r) where At 



is the typical spacing of the data points. 



We took the same set of 109 quasars that iKelly et al.l (120091 ) considered and fitted them 



with both approaches, finding results in reasonable agreement. We then took the parameters 
of these best fits and generated Monte Carlo realizations of the light curves (Equation ([3])) 
with the same time sampling and measurement uncertainties of the original light curves. We 
then re-fitted the Monte Carlo light curves using both methods and examined the differences 
between the input and output parameters. The results, as shown in FigureHU unambiguously 
show that the PRH method produces superior results, as would be expected from using all 
rather than only part of the correlation information in the data. 
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Fig. 1. — Examples of light curve models for four quasars. The 12- year-long light curves are 
from the OGLE-II and OGLE-III surveys (years 1997 - 2008, separated with the seasonal 
gaps). The top two light curves are sp ectroscopically confirmed quasars from Geha et al. 
(2003, top) and iDobrzvcki et al.l (120051) . The b ottom two light curves are mid-IR selected 
quasars from the iKozlowski fc Kochanekl (120091 ) sample (awaiting spectroscopic confirma- 
tion). The solid lines represent the best fit mean model light curves from the exponential 
covariance method (see appendix). The area between the dotted lines represent the la range 
of possible stochastic models. These "error snakes" bound the reconstructed light curve and 
are thinner than the data points because of the additional measurement error on the data 
(see the discussion of this point in Section [3]). We also give the best model parameter values 
along with the goodness of the fit defined by the x 2 P er degree of freedom. The average error 
bars on the model parameters are Alog 10 (r) = 0.37 and Alog 10 (<3-) = 0.044. 
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Fig. 2. — Distributions of the goodness of fit for the stochastic model (linear trend) are shown 
by the solid (dashed) line for the ~ 2, 700 OGLE-III light curves of the mid-IR selected 
quasars from iKozlowski &: Kochanekl (120091 ) . The improvement over fitting the linear trend 
from using the stochastic model is clearly visible. The dotted line is the expected distribution 
of x 2 /dof based on the number of degrees of freedom for each light curve. The differences 
between the stochastic model and the expected distribution are some combination of errors 
in the estimated photometric errors, small numbers of outliers that have not been eliminated 
from the light curves and any poorly modeled physics. 
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Fig. 3. — Fraction of variable QSO and YSO candidates as a function of magnitude for 
OGLE-III (solid line) and OGLE-II (dashed line). A source is counted as variable if it passes 
the variability criterion of Section [2] and has lnL^t > In L no i se + 2. 
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Fig. 4. — Monte Carlo simulations of completeness for source magnitudes of I = 16.5 (top 
left), 17.5 (top right), 18.5 (lower left) and 19.5 mag (lower right). The OGLE-II (solid line) 
and OLGE-III (dotted line) surveys have different cadences, durations and depths, leading 
to differences in their completeness limits. While not plotted, the distributions of variable 
sources at a given magnitude closely track these completeness limits in the sense that variable 
sources passing the selection criteria are not found in the regions where the completeness 
calculations say they should not be detected. The three lines show completeness levels of 
10%, 50%, and 90% (from left to right). The gray shaded region shows the region occupied 
by quasars (Cut 2 of Section [6]). 
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Fig. 5. — Locations of the ~ 1, 000 variable mid-IR selected quasar candidates from 



Kozlowski fc Kochanekl (120091 ) are shown in the CMD (top left), in aj-I space (top right), 



in a-T space (bottom left) and in (<jy/(j/)-(T/ space (bottom right). In the top left panel, the 
contours show the CMD of LMC stars, where the contours are for 1, 5, 10, and 20 stars per 
0.05 x 0.05 bin, counting from the outer contour. The RC is at (V — 1,1) = (1.05, 18.3) mag. 
In the lower le ft panel we show the line s of co nstant asymptotic variability a for 0.01, 0.1, and 
1.0 mag. The iKozlowski fc Kochanekl (120091 ) candidates are coded by color, where the high 
purity QSO-Aa objects are black, and the high contamination QSO-Ba and YSO-(AB)a are 
shown in red and green, respectively. The gray shaded regions indicate the regions occupied 
by quasars, and are defined by Cuts of Section [61 




Fig. 6. — Locations o f the ~ 1,000 variable mid-IR selected qu a sar ca ndidates from 



Kozlowski fc Kochanekl (120091 ) and ~ 100 quasars from iKelly et al.l (120091 ) are shown in 



<7-t space. The QSO-Aa objects are presented in the t op left pane l , QSO -Ba in the top right 
panel, YSO-(AB)a in the bottom right panel, and the lKelly et al.l ( 120091 ) quasars in the bot- 
tom left panel. The latter sample is split into Palomar Green (square), AGN Watch (star) 
and MACHO (triangle) quasars. We also show the lines of constant asymptotic variability 
a for 0.01, 0.1, and 1.0 mag (from left to right) and the region occupied by quasars (Cut 2 of 
Section [6j gray area). The thick solid lines in the bottom- left panel are for five AGN Watch 
quasars, and show the change in the model parameters after adding (subtracting) the host 
galaxy flux outside (inside) the spectroscopic aperture used for the light curve. Adding more 
host moves the source to longer time scales and lower amplitudes. 
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Fig. 7. — Our cleaned sample of ~ 10, 000 variable objects from the six central OGLE-II 
LMC fields (SCI - SC6) satisfying the basic variability criteria from Section [21 Panels are 
the same as in Figure [5] up to changes in the scales. The shaded regions show the quasar 
selection Cuts 2-4 from Section [6j 
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Fig. 8. — CMDs by variable type. In four panels, for clarity, we show distinct classes of 
variable objects. For a detailed description and discussion see text (Section EJ). Objects in 
the two masked areas in the left panels were moved from the bottom panel to the upper one. 
These objects were previously not cataloged by OGLE, but they are known variable classes. 
In the bottom-left panel, we show both periodic and non-periodic objects of unknown origin 
(or previously overlooked) and in the text we discuss the fraction that are quasars. The 
confirmed quasars and the OGLE-II quasar candidates have a similar distribution to the 
OGLE-III mid-IR selected quasar candidates (blue). For clarity, we do not include these 
smaller quasar samples here or in the subsequent figures. 
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Fig. 9. — Characteristic time as a function of modified amplitude (<t-t). The four panels 
show the same variability classes as in Figure The gray band is the quasar locus (Cut 2 
in Section [6]) used to separate quasars from the unknown non-periodic objects shown in the 
bottom-left panel. 
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Fig. 10. — Modified amplitude vs. magnitude (cx-I). The four panels show the same variability 
classes as in Figure El The gray area is the region defined by a high density of OGLE-III 
quasars, which separates them from other variability classes (Cut 3 in Section [6]). The 
heaviest contamination comes from the active giants and subgiants (AG/ASG, top left), and 
hot blue pulsators (SPB//3 Cep, top left) but these can be ruled out based on their periods 
(see Figure fT2l) . The gray wedge can in principle extend to arbitrarily faint magnitudes, but 
must be restricted based on the survey depth. For example, OGLE-II objects fainter than 
I > 19.5 mag are mostly false variables (below bottom horizontal dashed line). 
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Fig. 11. — Ratio of variability amplitudes in the V- and /-bands. It is clear that stellar 
objects generally vary more in the l/-band than in the /-band, while QSOs show smaller 
differences. The gray area is the region described by Cut 4 of Section [6j 
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period (days) 

Fig. 12. — Relation between Fourier periods P and time-scales r. These are not funda- 
mental and depend on the light curve sampling, particularly on the mean epoch spacing. 
Except for the bottom left panel, the panels show objects with a high probability of being 
periodic (log 10 (p per iodic) < —3) and non-periodic quasar candidates (blue dots, here selected 
with log 10 (p P eriodic) > — 3). Many periodic variable objects show a clear correlation of the 
characteristic time-scale r with period P. These are the Cepheids, OSARGs and long pe- 
riod variables. Quasars clearly show no such dependence. The bifurcation in period for 
periodic objects at r > 20 days (top panels) is probably due to the multi-mode pulsa- 
tions of OSARGs and LSPs. The solid line corresponds to a r-period relation described by 
log 10 (r) = 2 log 10 (P) — 1.8. The objects in the bottom left panel are previously unrecognized 
or overlooked non-periodic objects (black, log 10 (p per iodic) > — 3 or P > 200 days) and peri- 
odic variables (green, log 10 (p pe riodic) < —3 and P < 200 days). Spikes at 1, 1/2, 1/3 and 1/4 
days are due to aliasing problems in the periodograms that are not fully excluded by the 
narrow period notches we use to suppress most diurnal aliasing problems. 
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Fig. 13. — OGLE-II and OGLE-III light curve of the new variability-selected quasar can- 
didate LMC_SC5 253536. The variability behavior of the light curve closely resembles that 
of quasars, and it has met the basic variability criteria of Section [2] and all four cuts from 
Section [61 The solid line represents the best-fit model. As in Figure [Tj the area between the 
dotted lines represent the la range of possible stochastic models. We also give the process 
parameters and goodness of fit based on either the OGLE-II data or the combined data. 
The error bars on the model parameters are Alog 10 (r) = 0.49 and Alog 10 (cr) = 0.089. The 
vertical dashed line separates two phases of the OGLE survey. 
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Fig. 14. — Monte Carlo comparisons of p arameter est i matio n using the PRH and forecasting 
methods. For each of the 109 quasars in iKelly et al.l (120091 ) we generated four Monte Carlo 
realizations of each light curve using the either the best forecasting solution (open points) 
or PRH solution (dots) and then refit all eight of these light curve realizations using the 
forecasting (right) or PRH (left) method. The greater statistical power of the PRH approach 
is obvious. 
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Table 1. Sequential Cuts on Variable Objects. 



ij ui v ey 


\ /I ofriho/H 

IVldlCIlcU. 


Object 


Variable Not Noise 


Cut 1 






Pn+ ZL 




to OGLE 


Type 


(> 2a) 




Aperiodic 


<7-T 


a-I 








OGLE-III - Quasar Candidates 














/ > 16 mag 










OGLE-III 


2,156 


QSO-Aa 


721 


613 


598 


413 


401 


373 


OGLE-III 


215 


QSO-Ba 


75 


63 
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35 


35 


33 


OGLE-III 
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YSO-Aa 
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47 


28 
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14 


9 
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OGLE-III 


215 


QSO-Ba 


75 


63 


62 


35 


23 


23 


OGLE-III 
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OGLE-II 


/ only 


no masks 
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I only 
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13,658 
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V and / 
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/ only 
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11 
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V and I 
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Table 2. Individual Cuts on Variable Objects. 



Survey Matched Object Variable Not Noise Cut 1 Cut 2 Cut 3 Cut 4 
to OGLE Type (> 2a) Aperiodic o-t a-I oy/67 

OGLE-III - Quasar Candidates 



I > 16 mag 


OGLE-III 


2,156 


QSO-Aa 


721 


1,303 


2,131 


1,043 


1,560 


1,328 


OGLE-III 


215 


QSO-Ba 


75 


146 


231 


110 
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140 


OGLE-III 


198 


YSO-Aa 


147 


175 


169 


77 


84 


137 


OGLE-III 


24 


YSO-Ba 


15 


21 


23 


13 


14 


19 



OGLE-III - Quasar Candidates 
16 < / < 19.5 mag 



OGLE-III 


2,156 


QSO-Aa 
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1,303 


2,131 


1,043 
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1,328 
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63 
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21 


23 


13 


10 
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OGLE-II - Quasar Candidates 
16 < I < 19.5 mag 
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4 
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54 
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31 


48 


37 


15 


21 


19 
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3 


YSO-Ba 


1 


3 


3 





2 


1 



OGLE-II LMC Inner Fields - All Sources 
16 < / < 19.5 mag 

OGLE-II / only no masks 86,301 64,834 37,599 14,357 6,314 14,312 

OGLE-II / only with masks 13,658 10,406 3,461 1,980 800 5,433 

OGLE-II V and / no masks 71,247 53,942 30,798 12,182 5,333 14,312 

OGLE-II V and / with masks 12,729 9,973 3,219 1,853 706 5,433 



OGLE-II LMC Outer Fields - All Sources 
16 < / < 19.5 mag 



OGLE-II 


/ only 


no masks 


25,956 


19,251 


16,173 


2,508 


968 


2,056 


OGLE-II 


/ only 


with masks 


3,446 


1,272 
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OGLE-II 


V and / 


no masks 


16,404 


12,152 


9,820 


1,833 
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OGLE-II 


V and I 


with masks 


1,630 


443 


312 


98 


142 


70 



